function [mom, panel_out] = compute_moments(eq_SS, param, glob, options)

    %% markups
    eq_SS.mc     = param.omega*1./exp(glob.sf(:,2));
    eq_SS.markup = eq_SS.pf./eq_SS.mc;
    eq_SS.sales  = eq_SS.pf.*eq_SS.yf;
        
    eq_SS.costs  = eq_SS.l*param.omega;

    
    mom.mu_sw    = sum(eq_SS.sales.*eq_SS.L.*eq_SS.markup)/sum(eq_SS.sales.*eq_SS.L);
    mom.mu_cw    = sum(eq_SS.costs.*eq_SS.L.*eq_SS.markup)/sum(eq_SS.costs.*eq_SS.L);

    mom.mu_uw    = sum(eq_SS.L.*eq_SS.markup)/sum(eq_SS.L);
    
    q            = eq_SS.yf;
    elast        = param.sigma*q.^(-param.epsilon/param.sigma);
    mu_sigma     = elast./(elast-1);
    mom.cw_sigma = sum(mu_sigma.*eq_SS.costs.*eq_SS.L)/sum(eq_SS.costs.*eq_SS.L);
    
    mom.mu_disp  = sum((eq_SS.markup - mom.mu_uw).^2.*eq_SS.L);
    
    mom.mu_sales = sum(eq_SS.sales.*eq_SS.L);
    
    mom.sales_disp  = sum((eq_SS.sales - mom.mu_sales).^2.*eq_SS.L);
    
    param.w = 1; param.D = eq_SS.D; param.C = 1;
    
    acr             =  menufun('C',glob.sf,eq_SS.yf,param,glob,options); 
    
    mom.cost_size = sum(eq_SS.L.*acr)/sum(eq_SS.L.*eq_SS.sales);
    
    %%
     Lf = reshape(eq_SS.L, glob.Nsf,1);
     F = randsample(size(Lf,1), 1e6, true, Lf);

     prices  = eq_SS.pf(F);
     sales   = eq_SS.yf(F).*prices;
          
    prices  = eq_SS.pf(F);
    sales   = eq_SS.yf(F).*prices;

     
     prods   = glob.sf(F,2);

     markups = prices.*exp(prods);
     costs   = sales./prices./exp(prods);
 

    %% top 1% and 5% sales share
    mom.top10 = sum(maxk(sales, floor(.10*length(sales))))/sum(sales);
    mom.top05 = sum(maxk(sales, floor(.05*length(sales))))/sum(sales);
    mom.top01 = sum(maxk(sales, floor(.01*length(sales))))/sum(sales);

    mom.bottom87 = sum(mink(sales, floor(.877*length(sales))))/sum(sales); %15.4 percent of sales
    mom.frac_rel_sales_below_1 = sum(sales < mean(sales))/length(sales);

     %% simulate panel
    panel = ddpsimul(eq_SS.Q, F(1:100000),1);
  
%%
    
    prices  = eq_SS.pf(panel);

    prods   = glob.sf(:,2);
    prods   = prods(panel);
    labor   = eq_SS.yf(panel)./exp(prods);
    sales   = eq_SS.yf(panel).*prices;
        
    % I want to pick a sample that "looks like" Compustat.
    % If I pick big firms in first period, in 2nd period they're too small
    % because of mean reversion.
    % Instead, pick top firms among sum of total output in 1st and second
    % period
    

    sales   = sum(sales,2);
    labor   = sum(labor,2);
    [~,y] = sort(sales, 'descend');
    
    sample = labor(y);
    idx    = floor(length(sample)*.01); find(cumsum(sample) >= .3*sum(sample), 1);
    idx    = y(1:idx);
    
    sales   = eq_SS.yf(panel).*prices;
    labor   = eq_SS.yf(panel)./exp(prods);

    panel = panel(idx,:);
    prices  = prices(idx,:);
    prods   = prods(idx,:);
    sales   = sales(idx,:);
    markups = prices.*exp(prods);
    
    labor = labor(idx,:);
    
    mom.labor_churn = var(diff(log(labor),[],2));
    mom.sales_churn = var(diff(log(sales),[],2));
    mom.markups_churn = var(diff(log(markups),[],2));
    
    panel_out = [];
    panel_out.panel = panel;
    panel_out.labor = labor;
    panel_out.sales = sales;
    panel_out.markups = markups;
    panel_out.prices = prices;
    panel_out.prods = prods;
   
    
    %% run markups regression within firm
    %  also do labor regression
    
    %% compute slope coefficient
    lm = fitlm(log(sales(:,2)) - log(sales(:,1)),log(markups(:,2)) - log(markups(:,1)));
    mom.beta_mu = lm.Coefficients{2,1};
    
    lm = fitlm(log(sales(:,2)) - log(sales(:,1)),log(labor(:,2)) - log(labor(:,1)));
    mom.beta_l = lm.Coefficients{2,1};

    %%
    lm = fitlm(log(sales(:)),log(labor(:)));
    mom.beta_l_btw = lm.Coefficients{2,1};


    %% labor adjustment among large firms
    
    mom.large_jcjd10 = sum(diff(log(labor),[],2) > .1)/size(labor,1); % this moment is 28% in compustat (sample: post 2010)
    
    %% job creation/job destruction rates
    
    dL = (eq_SS.l - glob.sf(:,1))./glob.sf(:,1);
    %dL(glob.sf(:,1) == 0) = 0;
        
    mom.jc10 = sum(eq_SS.L.*(dL < .1).*(dL > 0))/sum(eq_SS.L);  % target: ~.05 quarterly, so ~.0125?
    mom.jd10 = sum(eq_SS.L.*(dL > -.1).*(dL < 0))/sum(eq_SS.L); % target: ~.05

    %% distribution of relative sales
    % whole economy statistics
         Lf = reshape(eq_SS.L, glob.Nsf,1);
     F = randsample(size(Lf,1), 10000, true, Lf);

     
    panel = ddpsimul(eq_SS.Q, F,2);
    
    labor = eq_SS.l(panel);
    
    sales = eq_SS.p(panel).*eq_SS.y(panel);
    
    d_labor = [labor(:,2) - labor(:,1), labor(:,3) - labor(:,2)];
    d_labor = d_labor./(.5*labor(:,1:2) + .5*labor(:,2:3));
    
    d_sales = [sales(:,2) - sales(:,1), sales(:,3) - sales(:,2)];
    d_sales = d_sales./(.5*sales(:,1:2) + .5*sales(:,2:3));
 
    
    mom.reallocation_rate = mean(mean(abs(d_labor)));
    
    labor_moments   = cov(d_labor); 
   	sales_moments   = cov(d_sales); 
    
    mom.demp_auto   = labor_moments(1,2)/labor_moments(1,1);
    mom.demp_var    = labor_moments(1,1);
    
    mom.dsales_auto = sales_moments(1,2)/labor_moments(1,1);
    mom.dsales_var  = sales_moments(1,1);
    
 %   mom.hiring_rate       = mean(d_labor);
    
    
    
end

